﻿using System;
using System.IO;
using System.IO.MemoryMappedFiles;
using System.Runtime.CompilerServices;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;

public struct PointStats
{
    public double Mean;
    public double Variance;
    public float Min;
    public float Max;
    public int Count;

    public override string ToString() =>
        $"Points: {Count}\n" +
        $"Mean: {Mean:F4}\n" +
        $"Variance: {Variance:F6}\n" +
        $"Min: {Min:F4}\n" +
        $"Max: {Max:F4}";
}

public static class PointStatsCalculator
{
    private const int ShuffleControlG4_YXWZ = 0b10110001; // 78: 交换高64位和低64位
    private const int ShuffleControl_ZWXY = 0b01001110; // 177: 交换每对64位内的两个32位元素

    private struct BatchStats
    {
        public float Min;
        public float Max;
        public double Sum;
        public double SumSquares;
        public int Count;
    }

    public static PointStats CalculateStats(string csvPath)
    {
        using var fileStream = new FileStream(csvPath, FileMode.Open, FileAccess.Read, FileShare.Read, 4096, FileOptions.SequentialScan);
        using var mmap = MemoryMappedFile.CreateFromFile(fileStream, null, 0, MemoryMappedFileAccess.Read, HandleInheritability.None, false);
        using var accessor = mmap.CreateViewAccessor(0, 0, MemoryMappedFileAccess.Read);

        return ParseAndCalculateStats(accessor);
    }

    private static unsafe PointStats ParseAndCalculateStats(MemoryMappedViewAccessor accessor)
    {
        const int batchSize = 8192;
        var valueBuffer = new float[batchSize];

        var globalStats = new BatchStats
        {
            Min = float.MaxValue,
            Max = float.MinValue,
            Sum = 0.0,
            SumSquares = 0.0,
            Count = 0
        };

        byte* ptr = null;
        accessor.SafeMemoryMappedViewHandle.AcquirePointer(ref ptr);

        try
        {
            long offset = 0;
            long size = accessor.Capacity;
            int bufferIndex = 0;
            int lineNumber = 0;
            int failedParseCount = 0;

            // 跳过UTF-8 BOM (如果存在)
            if (size >= 3 && ptr[0] == 0xEF && ptr[1] == 0xBB && ptr[2] == 0xBF)
            {
                offset += 3;
            }

            while (offset < size)
            {
                lineNumber++;
                long lineStartOffset = offset;

                // 解析一行中的所有数值
                int valuesParsed = ParseValuesInLine(ptr, ref offset, size, valueBuffer, bufferIndex, batchSize - bufferIndex);

                if (valuesParsed > 0)
                {
                    bufferIndex += valuesParsed;

                    // 当缓冲区满时，处理这一批数值
                    if (bufferIndex >= batchSize)
                    {
                        var batchStats = ProcessBatch(valueBuffer, bufferIndex);
                        MergeStats(ref globalStats, batchStats);
                        bufferIndex = 0;
                    }
                }
                else
                {
                    failedParseCount++;
                    SkipToLineEnd(ptr, ref offset, size);
                }
            }

            // 处理剩余的数值
            if (bufferIndex > 0)
            {
                var batchStats = ProcessBatch(valueBuffer, bufferIndex);
                MergeStats(ref globalStats, batchStats);
            }
        }
        finally
        {
            accessor.SafeMemoryMappedViewHandle.ReleasePointer();
        }

        return CalculateFinalStats(globalStats);
    }

    private static unsafe int ParseValuesInLine(byte* ptr, ref long offset, long size, float[] buffer, int startIndex, int maxValues)
    {
        int valuesParsed = 0;
        long lineStart = offset;

        // 跳过行首空白
        SkipWhitespace(ptr, ref offset, size);

        // 读取一行中的所有值
        while (offset < size && valuesParsed < maxValues)
        {
            // 尝试解析一个数值
            if (!TryParseFloat(ptr, ref offset, size, out float value))
            {
                break;
            }

            // 成功解析一个数值
            buffer[startIndex + valuesParsed] = value;
            valuesParsed++;

            // 跳过分隔符（准备解析下一个数值）
            SkipWhitespace(ptr, ref offset, size);
            if (offset < size && (ptr[offset] == ',' || ptr[offset] == ';' || ptr[offset] == '\t'))
            {
                offset++;
                SkipWhitespace(ptr, ref offset, size);
            }

            // 检查是否到达行尾
            if (offset >= size || ptr[offset] == '\n' || ptr[offset] == '\r')
            {
                break;
            }
        }

        // 跳过行尾
        SkipToLineEnd(ptr, ref offset, size);

        return valuesParsed;
    }

    private static BatchStats ProcessBatch(float[] values, int count)
    {
        // 使用SIMD优化处理，如果可用且数据量足够大
        if (Avx.IsSupported && count >= Vector256<float>.Count * 2)
        {
            return ProcessBatchAvx(values, count);
        }
        else if (Sse.IsSupported && count >= Vector128<float>.Count * 2)
        {
            return ProcessBatchSse(values, count);
        }
        else
        {
            return ProcessBatchScalar(values, count);
        }
    }

    private static unsafe BatchStats ProcessBatchAvx(float[] values, int count)
    {
        fixed (float* ptr = values)
        {
            int vectorSize = Vector256<float>.Count;
            int i = 0;

            // 初始化SIMD寄存器
            Vector256<float> minVec = Avx.LoadVector256(ptr);
            Vector256<float> maxVec = minVec;
            Vector256<float> sumVec = Vector256<float>.Zero;
            Vector256<float> sumSqVec = Vector256<float>.Zero;

            // 处理向量化部分
            for (i = vectorSize; i <= count - vectorSize; i += vectorSize)
            {
                Vector256<float> data = Avx.LoadVector256(ptr + i);

                minVec = Avx.Min(minVec, data);
                maxVec = Avx.Max(maxVec, data);
                sumVec = Avx.Add(sumVec, data);
                sumSqVec = Avx.Add(sumSqVec, Avx.Multiply(data, data));
            }

            // 水平归约
            float min = HorizontalMin(minVec);
            float max = HorizontalMax(maxVec);
            double sum = HorizontalSum(sumVec);
            double sumSq = HorizontalSum(sumSqVec);

            // 处理剩余部分
            for (; i < count; i++)
            {
                float value = values[i];
                min = Math.Min(min, value);
                max = Math.Max(max, value);
                sum += value;
                sumSq += value * value;
            }

            return new BatchStats
            {
                Min = min,
                Max = max,
                Sum = sum,
                SumSquares = sumSq,
                Count = count
            };
        }
    }

    private static unsafe BatchStats ProcessBatchSse(float[] values, int count)
    {
        fixed (float* ptr = values)
        {
            int vectorSize = Vector128<float>.Count;
            int i = 0;

            // 初始化SIMD寄存器
            Vector128<float> minVec = Sse.LoadVector128(ptr);
            Vector128<float> maxVec = minVec;
            Vector128<float> sumVec = Vector128<float>.Zero;
            Vector128<float> sumSqVec = Vector128<float>.Zero;

            // 处理向量化部分
            for (i = vectorSize; i <= count - vectorSize; i += vectorSize)
            {
                Vector128<float> data = Sse.LoadVector128(ptr + i);

                minVec = Sse.Min(minVec, data);
                maxVec = Sse.Max(maxVec, data);
                sumVec = Sse.Add(sumVec, data);
                sumSqVec = Sse.Add(sumSqVec, Sse.Multiply(data, data));
            }

            // 水平归约
            float min = HorizontalMin(minVec);
            float max = HorizontalMax(maxVec);
            double sum = HorizontalSum(sumVec);
            double sumSq = HorizontalSum(sumSqVec);

            // 处理剩余部分
            for (; i < count; i++)
            {
                float value = values[i];
                min = Math.Min(min, value);
                max = Math.Max(max, value);
                sum += value;
                sumSq += value * value;
            }

            return new BatchStats
            {
                Min = min,
                Max = max,
                Sum = sum,
                SumSquares = sumSq,
                Count = count
            };
        }
    }

    private static BatchStats ProcessBatchScalar(float[] values, int count)
    {
        float min = values[0];
        float max = values[0];
        double sum = values[0];
        double sumSq = values[0] * values[0];

        for (int i = 1; i < count; i++)
        {
            float value = values[i];
            min = Math.Min(min, value);
            max = Math.Max(max, value);
            sum += value;
            sumSq += value * value;
        }

        return new BatchStats
        {
            Min = min,
            Max = max,
            Sum = sum,
            SumSquares = sumSq,
            Count = count
        };
    }

    private static void MergeStats(ref BatchStats global, BatchStats batch)
    {
        if (global.Count == 0)
        {
            global = batch;
            return;
        }

        global.Min = Math.Min(global.Min, batch.Min);
        global.Max = Math.Max(global.Max, batch.Max);
        global.Sum += batch.Sum;
        global.SumSquares += batch.SumSquares;
        global.Count += batch.Count;
    }

    private static PointStats CalculateFinalStats(BatchStats stats)
    {
        if (stats.Count == 0)
            return new PointStats();

        double n = stats.Count;
        double mean = stats.Sum / n;
        double variance = stats.SumSquares / n - mean * mean;

        // 使用Welford算法修正数值精度
        // 对于大数据集，简单算法可能会有精度问题
        if (stats.Count > 1000)
        {
            variance = CorrectVariance(stats.Sum, stats.SumSquares, stats.Count);
        }

        return new PointStats
        {
            Mean = mean,
            Variance = Math.Max(0, variance),
            Min = stats.Min,
            Max = stats.Max,
            Count = stats.Count
        };
    }

    // 使用Welford算法计算更精确的方差
    private static double CorrectVariance(double sum, double sumSquares, int count)
    {
        if (count < 2) return 0;

        double mean = sum / count;
        double variance = (sumSquares - sum * mean) / (count - 1);
        return Math.Max(0, variance);
    }

    // SIMD水平归约辅助方法
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static float HorizontalMin(Vector256<float> vec)
    {
        Vector128<float> lower = vec.GetLower();
        Vector128<float> upper = vec.GetUpper();
        Vector128<float> min128 = Sse.Min(lower, upper);

        // 进一步归约到标量
        if (Sse.IsSupported)
        {
            Vector128<float> shuf = Sse.Shuffle(min128, min128, ShuffleControlG4_YXWZ);
            Vector128<float> min1 = Sse.Min(min128, shuf);
            shuf = Sse.Shuffle(min1, min1, ShuffleControl_ZWXY );
            Vector128<float> min2 = Sse.Min(min1, shuf);
            return min2.ToScalar();
        }

        // 回退到标量归约
        float min = float.MaxValue;
        for (int i = 0; i < Vector256<float>.Count; i++)
        {
            min = Math.Min(min, vec.GetElement(i));
        }
        return min;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static float HorizontalMax(Vector256<float> vec)
    {
        Vector128<float> lower = vec.GetLower();
        Vector128<float> upper = vec.GetUpper();
        Vector128<float> max128 = Sse.Max(lower, upper);

        // 进一步归约到标量
        if (Sse.IsSupported)
        {
            Vector128<float> shuf = Sse.Shuffle(max128, max128, ShuffleControlG4_YXWZ);
            Vector128<float> max1 = Sse.Max(max128, shuf);
            shuf = Sse.Shuffle(max1, max1, 0b01001110);
            Vector128<float> max2 = Sse.Max(max1, shuf);
            return max2.ToScalar();
        }

        // 回退到标量归约
        float max = float.MinValue;
        for (int i = 0; i < Vector256<float>.Count; i++)
        {
            max = Math.Max(max, vec.GetElement(i));
        }
        return max;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static double HorizontalSum(Vector256<float> vec)
    {
        Vector128<float> lower = vec.GetLower();
        Vector128<float> upper = vec.GetUpper();
        Vector128<float> sum128 = Sse.Add(lower, upper);

        // 进一步归约到标量
        if (Sse.IsSupported)
        {
            Vector128<float> shuf = Sse.Shuffle(sum128, sum128, ShuffleControlG4_YXWZ);
            Vector128<float> sum1 = Sse.Add(sum128, shuf);
            shuf = Sse.Shuffle(sum1, sum1, 0b01001110);
            Vector128<float> sum2 = Sse.Add(sum1, shuf);
            return sum2.ToScalar();
        }

        // 回退到标量归约
        double sum = 0;
        for (int i = 0; i < Vector256<float>.Count; i++)
        {
            sum += vec.GetElement(i);
        }
        return sum;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static float HorizontalMin(Vector128<float> vec)
    {
        if (Sse.IsSupported)
        {
            Vector128<float> shuf = Sse.Shuffle(vec, vec, ShuffleControlG4_YXWZ);
            Vector128<float> min1 = Sse.Min(vec, shuf);
            shuf = Sse.Shuffle(min1, min1, 0b01001110);
            Vector128<float> min2 = Sse.Min(min1, shuf);
            return min2.ToScalar();
        }

        // 回退到标量归约
        float min = float.MaxValue;
        for (int i = 0; i < Vector128<float>.Count; i++)
        {
            min = Math.Min(min, vec.GetElement(i));
        }
        return min;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static float HorizontalMax(Vector128<float> vec)
    {
        if (Sse.IsSupported)
        {
            Vector128<float> shuf = Sse.Shuffle(vec, vec, ShuffleControlG4_YXWZ);
            Vector128<float> max1 = Sse.Max(vec, shuf);
            shuf = Sse.Shuffle(max1, max1, 0b01001110);
            Vector128<float> max2 = Sse.Max(max1, shuf);
            return max2.ToScalar();
        }

        // 回退到标量归约
        float max = float.MinValue;
        for (int i = 0; i < Vector128<float>.Count; i++)
        {
            max = Math.Max(max, vec.GetElement(i));
        }
        return max;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static double HorizontalSum(Vector128<float> vec)
    {
        if (Sse.IsSupported)
        {
            Vector128<float> shuf = Sse.Shuffle(vec, vec, ShuffleControlG4_YXWZ);
            Vector128<float> sum1 = Sse.Add(vec, shuf);
            shuf = Sse.Shuffle(sum1, sum1, 0b01001110);
            Vector128<float> sum2 = Sse.Add(sum1, shuf);
            return sum2.ToScalar();
        }

        // 回退到标量归约
        double sum = 0;
        for (int i = 0; i < Vector128<float>.Count; i++)
        {
            sum += vec.GetElement(i);
        }
        return sum;
    }

    #region CSV解析工具方法

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static unsafe void SkipWhitespace(byte* ptr, ref long offset, long size)
    {
        while (offset < size)
        {
            byte c = ptr[offset];
            if (c != ' ' && c != '\t' && c != '\r' && c != '\n')
                break;
            offset++;
        }
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static unsafe void SkipToLineEnd(byte* ptr, ref long offset, long size)
    {
        while (offset < size && ptr[offset] != '\n' && ptr[offset] != '\r')
            offset++;
        SkipLineEnding(ptr, ref offset, size);
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static unsafe void SkipLineEnding(byte* ptr, ref long offset, long size)
    {
        if (offset < size && ptr[offset] == '\r')
            offset++;
        if (offset < size && ptr[offset] == '\n')
            offset++;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static unsafe bool TryParseFloat(byte* ptr, ref long offset, long size, out float result)
    {
        result = 0;
        if (offset >= size) return false;

        // 跳过前导空格
        while (offset < size && (ptr[offset] == ' ' || ptr[offset] == '\t'))
            offset++;

        if (offset >= size) return false;

        bool negative = false;
        if (ptr[offset] == '-')
        {
            negative = true;
            offset++;
        }
        else if (ptr[offset] == '+')
        {
            offset++;
        }

        if (offset >= size) return false;

        bool hasDigits = false;
        bool hasDecimal = false;
        float decimalMultiplier = 0.1f;
        long integerPart = 0;

        // 解析整数部分
        while (offset < size)
        {
            byte c = ptr[offset];
            if (c >= '0' && c <= '9')
            {
                hasDigits = true;
                integerPart = integerPart * 10 + (c - '0');
                offset++;
            }
            else
            {
                break;
            }
        }

        result = integerPart;

        // 解析小数部分
        if (offset < size && ptr[offset] == '.')
        {
            hasDecimal = true;
            offset++;

            while (offset < size)
            {
                byte c = ptr[offset];
                if (c >= '0' && c <= '9')
                {
                    result += (c - '0') * decimalMultiplier;
                    decimalMultiplier *= 0.1f;
                    offset++;
                }
                else
                {
                    break;
                }
            }
        }

        // 解析指数部分
        if (offset < size && (ptr[offset] == 'e' || ptr[offset] == 'E'))
        {
            offset++;
            bool expNegative = false;

            if (offset < size && ptr[offset] == '-')
            {
                expNegative = true;
                offset++;
            }
            else if (offset < size && ptr[offset] == '+')
            {
                offset++;
            }

            int exponent = 0;
            while (offset < size && ptr[offset] >= '0' && ptr[offset] <= '9')
            {
                exponent = exponent * 10 + (ptr[offset] - '0');
                offset++;
            }

            if (exponent > 0)
            {
                float power = MathF.Pow(10, exponent);
                if (expNegative)
                    result /= power;
                else
                    result *= power;
            }
        }

        if (!hasDigits) return false;
        if (negative) result = -result;
        return true;
    }

    #endregion
}

class Program
{
    static void Main(string[] args)
    {
        string csvPath = "d:\\download\\sample-data.csv";
        if (!File.Exists(csvPath))
        {
            Console.WriteLine($"Error: File '{csvPath}' not found");
            return;
        }

        try
        {
            Console.WriteLine($"Processing file: {csvPath}");
            Console.WriteLine($"File size: {new FileInfo(csvPath).Length / 1024.0 / 1024.0:F2} MB");

            var stopwatch = System.Diagnostics.Stopwatch.StartNew();
            var stats = PointStatsCalculator.CalculateStats(csvPath);
            stopwatch.Stop();

            Console.WriteLine("\nStatistics:");
            Console.WriteLine(stats.ToString());
            Console.WriteLine($"\nProcessing time: {stopwatch.ElapsedMilliseconds}ms");

            if (stats.Count > 0)
            {
                double throughput = stats.Count / stopwatch.Elapsed.TotalSeconds;
                Console.WriteLine($"Throughput: {throughput:F0} points/second");

                // 检查SIMD支持
                Console.WriteLine($"\nHardware support:");
                Console.WriteLine($"AVX: {Avx.IsSupported}");
                Console.WriteLine($"SSE: {Sse.IsSupported}");
                Console.WriteLine($"Vector256: {Vector256<float>.IsSupported}");
                Console.WriteLine($"Vector128: {Vector128<float>.IsSupported}");
            }
        }
        catch (Exception ex)
        {
            Console.WriteLine($"Error processing file: {ex.Message}");
            Console.WriteLine($"Stack trace: {ex.StackTrace}");
        }
    }
}